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Abstract 


The  time  of  flight  of  a  two-body  orbit  may  be  determined  by  integrating  the  radial 
velocity  equation  for  a  conic  section.  The  resulting  expression  is  sometimes  called  Lam¬ 
bert’s  Time  Function,  which  depends  on  the  gravitational  constant,  two  position  vectors, 
and  the  semi-major  axis  of  the  conic  flight  path.  For  mission  planning  purposes,  it  is  often 
desirable  to  know  the  semi-major  axis  as  a  function  of  time,  rather  than  the  reverse.  Nor¬ 
mally,  a  root  finding  technique  such  as  Newton- Raphson  is  employed  to  find  the  value  of  a 
characteristic  orbital  parameter  which  matches  a  given  time  of  flight.  Alternatively,  Lam¬ 
bert’s  Time  vunction  may  be  expanded  as  a  power  series  involving  the  inverse  semi-major 
axis.  The  expression  for  semi-major  axis  is  then  determined  through  series  reversion  and 
inversion  of  the  resulting  series.  A  simplified  method  of  obtaining  the  series  coefficients  is 


given,  as  well  as  a  numerical  study  of  convergence  properties. 


SERIES  REVERSION/INVERSION  OF 
LAMBERT’S  TIME  FUNCTION 

I.  Introduction 

The  two-position,  time  of  flight  problem  in  orbital  mechanics  has  historically  been 
the  subject  of  many  investigations.  Although  a  closed  form  solution  is  available  through 
integration  of  the  polar  velocity  equation  for  a  conic  section,  the  resulting  expression  is 
transcendental  in  the  orbital  parameters,  making  it  difficult  to  solve  for  them.  In  order  to 
match  the  orbital  parameters  to  a  particular  flight  time,  root  finding  techniques  are  nor¬ 
mally  used.  Many  different  formulations  of  the  problem  have  been  developed  to  minimize 
the  convergence  time  of  these  root  finding  methods,  as  well  as  to  generalize  the  problem 
to  avoid  case  dependent  equations.  In  each  of  these  methods,  an  initial  value  is  required. 
Depending  on  the  method  used,  convergence  may  not  be  achieved  at  all  if  the  initial  value 
is  too  different  from  the  correct  value.  In  the  classical  Gauss  method  (4:  p.  188-197  ), 
for  example,  the  algorithm  fails  to  converge  for  transfer  angles  larger  than  roughly  70°. 
Although  other  methods  have  been  developed  to  converge  for  larger  angles,  they  may  be 
deficient  for  small  angles.  It  would  be  desirable,  therefore,  to  develop  a  solution  that 
does  not  require  an  initial  value  and  is  not  dependent  on  transfer  angle  for  convergence 
properties. 

The  method  presented  is  a  solution  of  the  closed  form  time  of  flight  equation  for 
semi-major  axis  as  a  function  of  time.  It  will  be  shown  that  this  function  may  be  used  to 
directly  calculate  the  semi-major  axis  of  an  orbit  without  the  utilization  of  a  root  finding 
technique. 
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II.  Analytical  Development 


The  Two-Body  Time  of  Flight  Problem 

The  two-body  time  of  flight  problem  may  be  stated  as  follows:  Given  two  position 
vectors  from  the  gravitational  body  to  the  orbiting  body  and  :■  time  of  flight  between  the 
positions,  determine  the  semi-major  axis  of  the  conic  trajectory  followed  by  the  body. 


Figure  1.  Geometry  of  two-body  problem 


In  the  problem  statement  above  it  is  assumed  that  the  Newtonian  gravitational  con¬ 
stant  is  known,  and  both  bodies  behave  as  point  masses. 
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Lambert’s  Time  Function 


The  purpose  of  this  section  is  to  derive  Lambert’s  Time  Function  in  a  closed,  tran¬ 
scendental  form.  Once  known,  it  can  then  be  expressed  as  a  power  series  allowing  for  later 
reversion  and  inversion  to  solve  for  semi-major  axis. 

If  the  position  of  the  orbiting  body  is  expressed  in  polar  coordinates,  the  energy 
integral  may  be  written: 


(5)* -(H) 

where  r  is  the  radial  distance  between  the  two  bodies,  n  is  the  gravitational  constant, 
a  is  the  semi-major  axis  of  the  conic  section,  and  t  is  the  time.  Equation  (1)  may  be 
rearranged  as  follows: 


dt  = 


r  dr 


</jiy/2r-(r*/a) 
The  integral  takes  the  form  (3:  p.  70-75  ): 


(2) 


-f 

Ja-c 


rdr 


-C  y/JI^r^Jr^Ja) 

where  s  =  (ri  +  r2  +  c)/2,  and  c  =  chord.  (  see  figure  1.) 

The  integration  may  be  simplified  by  introducing  a  change  of  variable: 


(3) 


r  =  a(l  -  cos<f>) 

To  find  the  new  limits  of  integration: 

ri  =  s  —  c  =  a(l  -  cos  <j> ) 


(4) 

(5) 


<t>  =  cos  1  (1  -  - — -) 
a 


(6) 
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Similarly,  the  upper  limit  of  integration  becomes: 


t  =  ^-[(a- sin a)-(/? -sin/?)]  (12) 

This  is  Lambert’s  Time  Function  for  an  elliptic  trajectory  with  a  transfer  angle  less  than  7r, 
and  a  flight  time  less  than  the  minimum  energy  transfer  time.  A  minimum  energy  transfer 
is  an  elliptical  transfer  where  the  semi-major  axis  is  equal  to  one  half  of  the  semi-perimeter, 
a  =  s/2. 
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Series  Representation  of  Time  of  Flight 

In  order  to  accomplish  series  reversion  and  inversion,  the  time  of  flight  must  first  be 
expressed  as  a  single  power  series.  Substitution  of  the  definitions  of  a  and  / 3  yields: 


Substituting  into  equation  (13)  produces: 


t 


(i  Ll.±)_x[±F(zl 
\2’2’2’2  a)  V  2a  \  2 


5 


l/Sfe)"  Kf  H;  i)  - r  (T";i)] - 


2/i 


(14) 


The  hypergeometric  functions  may  be  expressed  as  infinite  series,  allowing  manipula¬ 
tion  of  individual  terms  to  form  a  single,  combined  series.  Representing  the  hypergeometric 
functions  with  Pochhammer  notation  yields: 


The  first  terms  in  both  series  of  equation  (15)  are  zero.  Changing  the  indicies  by  m  =  n- 1, 
n  =  m  +  1  produces: 


t  = 


0) 


.  o3  oo 

ll_  y* 

2M  ~ 
r  m=0 


(«  -  03  v 


2 

3  V  2m 

*  ^  m=0 


m  +  1 


m+1  ((rn  +  l)2  -  1/4) 
m  -(-  1 


(0 


m+1  ((m+1)2  -  1/4) 

(I). 


l  .(±y_ 

(m  +  1)!  \2aJ 

-  I—  flzfY 

(m  +  1)!  \  2a  ) 


(4/3)(m  +  l)2  -  1/3 


(m)!  \2a/ 


2  /(,  -  c)3  - 

3  V  2m  ^ 


m=0 


1  ' 

fa> 

v2y 

1 

_ 

i 

(3~C\ 

[(4/3)(m+l)2-l/3_ 

(m)! 

\  2  a  ) 

(16) 
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In  order  to  produce  the  hypergeometric  parameters  (1:  p.  276-277),  one  must  form 
the  ratio  u^‘  where  um  is  the  mth  term  of  the  infinite  series.  This  gives  the  hypergeometric 
parameters  (5,  §,  §)  by  inspection: 


(um+i/um)  = 


'(m+  l/2)(m  4-  3/2) 
(m  +  5/2) 


1 

m  +  1 


Substitution  into  equation  (16)  yields: 


(17) 


The  first  term  may  be  recognized  as  the  Parabolic  Transfer  Time,  tp  ,  for  the  given 
problem  geometry.  (  see  figure  1  )  The  Parabolic  Transfer  Time  is  the  flight  time  for  an 
object  traveling  on  a  parabolic  trajectory  between  two  position  vectors,  where  the  origin 
of  the  position  vectors  is  the  focus  of  the  parabola.  The  flight  time  of  an  object  traveling 
on  a  hyperbolic  trajectory  will  be  less  than  the  Parabolic  Transfer  Time,  and  the  flight 
time  on  an  elliptic  trajectory  will  be  greater.  The  semi-major  axis  of  a  parabola  is  infinite 
by  definition,  and  is  the  limiting  case  between  hyperbolas  and  ellipses.  Let  the  quantity 
^  -  l)  =  T,  a  non-dimensional  time  parameter. 
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(19) 


The  gravitational  constant  n  is  included  only  in  the  non-dimensional  time  parameter  T. 
Therefore,  the  series  coefficients  are  functions  of  geometric  constants. 


It  may  be  shown  that  the  time  equation  for  the  hyperbolic  case  differs  only  in  the 
sign  of  the  argument,  i.  e.,  (ff),  and  the  series  coefficients  are  identical. 
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Reversion/ Inversion  of  Series 


The  purpose  of  this  section  is  to  solve  equation  (19)  for  the  semi-major  axis,  a.  To 
this  end,  a  series  reversion  will  be  accomplished  followed  by  a  binomial  expansion  of  its 
inverse.  In  general,  equation  (19)  may  be  expressed  as  follows: 

T = M  (£) +  (£)  +/l3(i)  +-  (20) 

In  a  series  reversion,  the  expression  is  rewritten  (2:  p.  316-317): 

(i)-4T+4r*+4r*+_  (21) 

In  order  to  determine  the  semi-major  axis,  the  inverse  must  be  found: 

(7)  =  (a[t  +  a'2t2  +  a'3t3  +  ...y1 

=  £(-1  )n(A'1T)(-1-n)(A;r2  +  a'3t3  +  ...)n 

n=0 

=  (A'/ry1  -  (A[Ty\A'2T2  +  a3t3  + ...) 

+{A'lTy3(Al2T2  +  A3T3  +  ...)2  -  ...  (22) 

Or  more  generally: 

=  BxT~l  +  B2T°  +  B3Tx  +  B4T 2  +  ...  (23) 

by  equating  powers  of  T. 
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Determination  of  Series  Coefficients 


Now  that  the  semi-major  axis  has  been  expressed  as  a  power  series  in  T,  a  method 
must  be  developed  to  determine  the  unknown  coefficients  B,  in  terms  of  the  known  coeffi¬ 
cients  A,.  Equation  (23)  is  repeated  here  for  convenience: 

(y)  =  BXT-1  +  B2T°  +  B3T1  +  B4T 2  +  ...  (24) 


Multiply  by  (^)  T  and  substitute  for  T: 

Al  (3j)  +  a2  (^)2  +  A3  (^)3  +  ...  = 

(^)  +  -®2  (^)  (>li  (5;)  +  A2  (^)2  +  A3  (^)3  +  ...}  +  ... 

Taking  a  derivative  with  respect  to  (^)  yields: 

Ai  +  2A2  (£)  +  3A3  (^)2  +  ...  = 

ft  +  ft  {a.  (4)  +  ft  (4)!  +  ft  (4)3  +  ■■•}  (26) 

+ft  (4)  {ft  +  2ft  (4)  +  3ft  (4)!  + ...}  + ... 


Evaluate  at  (^)  =  0: 

Ai  =  B\  (27) 

In  an  analogous  fashion,  successive  derivatives  evaluated  at  zero  produce  the  follow¬ 

ing: 

A2  =  A\B2  (28) 

A3  =  A2B2  +  A2B3  (29) 

The  expressions  become  increasingly  complex  with  each  derivative. 

In  general,  a  matrix  equation  may  be  formed  showing  the  relationships  between  A, 
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and  B{: 


a2 

O 

O 

O 

O 

H 

1 _ 

Bi  ' 

A3 

A2  A\  0  0  0 

S3 

a4 

►  = 

A3  2  A1A2  A\  0  0 

< 

s4 

■ 

.  .  0 

• 

Ai 

.  .  •  A\ 

Bi 

Solving  for  the  f?,: 


Q~lS 


(30) 


B2 

S3 

-a2/a\ 

S4 

>  = 

(2 A\  -  AiA3)/A 

Bi 

To  summarize: 


0 

0 

0 

0 

’  > 

a2 

Ar2 

0 

0 

0 

a3 

2A2/A? 

Aj-3 

0 

0 

- 

a4 

• 

• 

. 

0 

■ 

• 

• 

• 

A,- 

(31) 


Bi  = 

52  =  A2/A1 

53  =  (AxAz-AD/Al 

54  —  ( A2A4  —  3Ai  A2A3  +  2A%)/A\ 


(32) 

(33) 

(34) 

(35) 


It  may  be  seen  that  the  elements  of  the  matrix  in  equation  (31)  follow  a  pattern.  This 
allows  the  matrix  to  be  formed  without  the  necessity  of  taking  numerous  derivatives  as  was 
shown  previously.  As  an  aid  to  finding  this  pattern,  the  computer  program  MACSYMA,  a 
trademark  of  Symbolics,  Inc.  (Project  MAC’s  SYmbolic  MAnipulation  System)  was  used 
to  generate  numerous  terms  for  inspection.  The  infinite  series  were  represented  by  tenth 
order  polynomials  and  then  expanded  by  MACSYMA,  which  also  evaluated  ten  derivatives 
in  the  manner  of  equation  (26).  This  produced  ten  relationships  between  the  Ai  and  the 
Bi  which  were  arranged  in  matrix  form.  Finally,  MACSYMA  was  used  to  symbolically 
invert  a  subset  of  the  matrix  allowing  for  the  pattern  to  be  recognized  through  inspection, 
trial  and  error.  It  should  be  emphasized  that  a  symbolic  manipulation  program  such  as 
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MACSYMA  may  prove  to  be  an  invaluable  tool  in  this  type  of  analysis. 

To  form  the  matrix  in  equation  (31),  element  q\\  is  assumed  to  be  AJ-1.  From  this, 
all  remaining  elements  may  be  formed.  The  process  is  illustrated  by  example: 

To  form  element  933,  take  the  following  dot  product: 


933  =  {921,911}  •  {912, 922}  (36) 

=  {-A2M?,Arx}.{o,Ar2}  =  V  (37) 

To  form  element  932,  take  the  following  dot  product: 

932  =  {921  >9n}  •  {911  >  921}  (38) 

=  {-A2/A?,Ar1}.{Aj-1,-A2/A?}  =  — 2A2/Aj  (39) 


To  form  element  931  : 

931  =  -~r~  {932,933,934>—}  •  {^2,  A3,  A4, ...} 

Ai 

=  ~  {-2At/A*  +  A3^}  =  {2A2  -  MAz){A\ 


(40) 

(41) 


All  of  the  elements  of  the  third  row  are  now  determined,  and  the  fourth  row  may  be 
found  beginning  with  and  ending  with  941.  The  first  element  of  each  new  row  must 
be  calculated  last,  since  all  of  the  remaining  row  elements  are  required.  To  generalize  the 
procedures,  the  following  steps  should  be  performed: 

1-  9n  =  V 

2-  9u  =  (9(«-i,i)»9(«-2,i)>—9ii)  •  (9(ij'-i)»9(2,i-i)»  — 9(«-i,j-i))>  2  <  j  <=  i 

3.  qn  =  (-Af1){(g,-2,gi3,...9iii)«(-A71)(A2,A3,...A,)} 


4.  B{  =  (9ii.9i2>  "m9«»)  •  (A2,  A3, ...,  A,-+i) 
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Repeat  steps  2,  3,  and  4  for  each  additional  term  desired.  (  see  Appendix  C,  "Form 
Q  Matrix”  )  Because  the  A,  are  known  from  equation  (19)  they  may  be  numerically 
determined  beforehand  and  used  to  calculate  the  £,  coefficients. 

In  the  figure  below,  rx  =  r2  =  k  =  1,  and  6  =  180°.  Recall  T  =  0  corresponds  to  a 
parabolic  transfer. 


Figure  2.  Semi-major  axis  vs.  T  for  general  180°  case. 
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III.  Numerical  Investigation  and  Results 


General  Case 

The  general  case  considered  was  composed  of  the  following  orbital  geometry: 

•  magnitudes  of  position  vectors:  r\  =  r^  —  1 

•  gravitational  constant:  fi  =  1 

•  range  of  transfer  angle:  0°  <  9  <  360°,  15°  increments 

•  range  of  flight  time:  0  <  t  <  tme 

The  following  data  consists  of  values  of  semi-major  axis  calculated  for  each  15°  of 
transfer  angle  using  partial  sums  of  the  series.  Values  are  tabulated  for  terms  i  =  1,2,3 
and  i  =  21,22,23.  The  last  column  is  the  value  of  T  obtained  by  using  the  first  23  terms 
to  find  the  semi-major  axis,  then  by  employing  Lambert’s  Time  Function.  These  values  of 
T  may  be  used  as  an  accuracy  check  against  the  desired  value  of  T  shown  in  the  caption 
of  each  table.  Taking  an  example  row  of  data,  for  T  =  -0.9  and  9  —  30°,  the  semi-major 
axis  value  using  one  term  in  the  series  is  -0.28090,  for  two  terms  is  -0.10653,  and  for  terms 
21,  22,  and  23  remains  a  constant  -0.00504.  At  this  point  the  semi-major  axis  seems  to 
have  converged  to  the  correct  value,  since  application  of  Lambert’s  Time  Function  yields 
T  =  -0.9  as  the  last  entry  in  the  data  row. 

The  accuracy  of  the  series  is  plotted  against  T  and  i  in  Appendix  B  for  transfer 
angles  of  90°,  180°,  270°,  and  360°.  The  vertical  axis  represents  the  number  of  significant 
figures  in  the  semi-major  axis  value.  The  four  plots  have  a  very  similar  appearance,  which 
would  seem  to  indicate  a  general  insensitivity  of  accuracy  with  respect  to  the  transfer 
angle.  Plots  were  constructed  for  every  10°,  but  because  of  the  similarities,  only  four  plots 
are  presented. 
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Table  1.  Values  of  semi-major  axis  for  the  general  case,  T  =  -0.9 


e 

1 

2 

3 

15.0 

-0.27857 

0.09962 

0.03942 

30.0 

-0.28090 

0.10653 

0.03527 

45.0 

-0.28468 

0.11718 

0.02979 

60.0 

-0.28972 

0.13043 

0.02450 

75.0 

-0.29578 

0.14494 

0.02058 

90.0 

-0.30255 

0.15939 

0.01853 

105.0 

-0.30966 

0.17263 

0.01819 

120.0 

-0.31667 

0.18380 

0.01894 

135.0 

-0.32309 

0.19239 

0.02007 

150.0 

-0.32839 

0.19822 

0.02102 

165.0 

-0.33200 

0.20142 

0.02155 

180.0 

-0.33333 

0.20238 

0.02168 

195.0 

-0.33181 

0.20161 

0.02163 

210.0 

-0.32692 

0.19965 

0.02160 

225.0 

-0.31824 

0.19693 

0.02161 

240.0 

-0.30556 

0.19363 

0.02130 

255.0 

-0.28894 

0.18952 

0.01986 

270.0 

-0.26888 

0.18382 

0.01617 

285.0 

-0.24640 

0.17528 

0.00944 

300.0 

-0.22310 

0.16254 

0.00059 

315.0 

-0.20114 

0.14513 

-0.00620 

330.0 

-0.18296 

0.12509 

-0.00486 

345.0 

-0.17090 

0.10804 

0.00470 

360.0 

-0.16667 

0.10119 

0.01084 

21 

22 

23 

T 

-0.00504 

-0.00504 

-0.00504 

-0.9000 

-0.00502 

-0.00502 

-0.00502 

-0.9000 

-0.00499 

-0.00499 

-0.00499 

-0.9000 

-0.00494 

-0.00494 

-0.00494 

-0.9000 

-0.00489 

-0.00489 

-0.00489 

-0.9000 

-0.00483 

-0.00483 

-0.00483 

-0.9000 

-0.00476 

-0.00476 

-0.00476 

-0.9000 

-0.00470 

-0.00470 

-0.00469 

-0.9000 

-0.00464 

-0.00464 

-0.00464 

-0.8999 

-0.00458 

-0.00458 

-0.00458 

-0.9001 

-0.00456 

-0.00456 

-0.00456 

-0.9000 

-0.00456 

-0.00456 

-0.00456 

-0.9000 

-0.0P':. 

-0.00452 

-0.00453 

-0.9000 

-0.00439 

-0.00439 

-0.00439 

-0.9000 

-0.00415 

-0.00415 

-0.00415 

-0.9002 

-0.00393 

-0.00391 

-0.00391 

-0.8997 

-0.00356 

-0.00354 

-0.00354 

-0.9005 

-0.00335 

-0.00335 

-0.00332 

-0.8995 

-0.00287 

-0.00301 

-0.00311 

-0.8982 

-0.00249 

-0.00261 

-0.00288 

-0.8977 

-0.00221 

-0.00250 

-0.00283 

-0.8948 

-0.00233 

-0.00260 

-0.00232 

-0.9017 

-0.00233 

-0.00229 

-0.00233 

-0.8995 

-0.00228 

-0.00228 

-0.00228 

-0.9000 
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Table  2.  Values  of  semi-major  axis  for  the  general  case,  T  =  -0.4 


6 

1 

2 

3 

15.0 

-0.62678 

-0.24859 

-0.27535 

30.0 

-0.63203 

-0.24460 

-0.27627 

45.0 

-0.64053 

-0.23867 

-0.27751 

60.0 

-0.65187 

-0.23172 

-0.27880 

75.0 

-0.66550 

-0.22479 

-0.28006 

90.0 

-0.68074 

-0.21880 

-0.28140 

105.0 

-0.69673 

-0.21444 

-0.28308 

120.0 

-0.71250 

-0.21203 

-0.28530 

135.0 

-0.72695 

-0.21147 

-0.28806 

150.0 

-0.73888 

-0.21227 

-0.29102 

165.0 

-0.74700 

-0.21358 

-0.29352 

180.0 

-0.75000 

-0.21429 

-0.29460 

195.0 

-0.74658 

-0.21316 

-0.29315 

210.0 

-0.73557 

-0.20901 

-0.28814 

225.0 

-0.71604 

-0.20087 

-0.27879 

240.0 

-0.68750 

-0.18831 

-0.26490 

255.0 

-0.65011 

-0.17165 

-0.24706 

270.0 

-0.60498 

-0.15228 

-0.22679 

285.0 

-0.55439 

-0.13272 

-0.20642 

300.0 

-0.50198 

-0.11634 

-0.18831 

315.0 

-0.45256 

-0.10629 

-0.17355 

330.0 

-0.41165 

-0.10360 

-0.16136 

345.0 

-0.38452 

-0.10558 

-0.15151 

360.0 

-0.37500 

-0.10714 

-0.14730 

21 

22 

23 

T 

-0.28128 

-0.28128 

-0.28128 

-0.4000 

-0.28140 

-0.28140 

-0.28140 

-0.4000 

-0.28167 

-0.28167 

-0.28167 

-0.4000 

-0.28217 

-0.28217 

-0.28217 

-0.4000 

-0.28303 

-0.28303 

-0.28303 

-0.4000 

-0.28437 

-0.28437 

-0.28437 

-0.4000 

-0.28629 

-0.28629 

-0.28629 

-0.4000 

-0.28883 

-0.28883 

-0.28883 

-0.4000 

-0.29184 

-0.29184 

-0.29184 

-0.4000 

-0.29495 

-0.29495 

-0.29495 

-0.4000 

-0.29752 

-0.29752 

-0.29752 

-0.4000 

-0.29862 

-0.29862 

-0.29862 

-0.4000 

-0.29716 

-0.29716 

-0.29716 

-0.4000 

-0.29211 

-0.29211 

-0.29211 

-0.4000 

-0.28276 

-0.28276 

-0.28276 

-0.4000 

-0.26889 

-0.26889 

-0.26889 

-0.4000 

-0.25104 

-0.25104 

-0.25104 

-0.4000 

-0.23052 

-0.23052 

-0.23052 

-0.4000 

-0.20921 

-0.20921 

-0.20921 

-0.4000 

-0.18920 

-0.18920 

-0.18920 

-0.4000 

-0.17228 

-0.17228 

-0.17228 

-0.4000 

-0.15964 

-0.15964 

-0.15964 

-0.4000 

-0.15190 

-0.15190 

-0.15190 

-0.4000 

-0.14931 

-0.14931 

-0.14931 

-0.4000 
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Table  3.  Values  of  semi-major  axis  for  the  general  case,  T  —  0.1 


6  1 
15.0  2.50711 
30.0  2.52813 
45.0  2.56211 
60.0  2.60747 
75.0  2.66201 
90.0  2.72295 
105.0  2.78692 
120.0  2.85000 
135.0  2.90781 
150.0  2.95552 
165.0  2.98801 
180.0  3.00000 
195.0  2.98633 
210.0  2.94229 
225.0  2.86418 
240.0  2.75000 
255.0 
270.0 
285.0 
300.0 
315.0 
330.0 
345.0 
360.0 


2 

3 

2.88530 

2.89199 

2.91557 

2.92349 

2.96397 

2.97368 

3.02762 

3.03939 

3.10273 

3.11655 

3.18489 

3.20054 

3.26920 

3.28636 

3.35047 

3.36879 

3.42328 

3.44243 

3.48212 

3.50181 

3.52143 

3.54141 

3.53571 

3.55579 

3.51975 

3.53975 

3.46886 

3.48864 

3.37935 

3.39883 

3.24919 

3.26834 

3.07889 

3.09774 

2.87260 

2.89123 

2.63925 

2.65767 

2.39355 

2.41155 

2.15650 

2.17332 

1.95465 

1.96909 

1.81700 

1.82848 

1.76786 

1.77790 

21 

22 

2.89171 

2.89171 

2.92326 

2.92326 

2.97351 

2.97351 

3.03924 

3.03924 

3.11641 

3.11641 

3.20039 

3.20039 

3.28619 

3.28619 

3.36860 

3.36860 

3.44223 

3.44223 

3.50160 

3.50160 

3.54120 

3.54120 

3.55558 

3.55558 

3.53954 

3.53954 

3.48843 

3.48843 

3.39862 

3.39862 

3.26813 

3.26813 

3.09753 

3.09753 

2.89102 

2.89102 

2.65749 

2.65749 

2.41145 

2.41145 

2.17337 

2.17337 

1.96925 

1.96925 

1.82852 

1.82852 

1.77779 

1.77779 

23 

T 

2.89171 

0.1000 

2.92326 

0.1000 

2.97351 

0.1000 

3.03924 

0.1000 

3.11641 

0.1000 

3.20039 

0.1000 

3.28619 

0.1000 

3.36860 

0.1000 

3.44223 

0.1000 

3.50160 

0.1000 

3.54120 

0.1000 

3.55558 

0.1000 

3.53954 

0.1000 

3.48843 

0.1000 

3.39862 

0.1000 

3.26813 

0.1000 

3.09753 

0.1000 

2.89102 

0.1000 

2.65749 

0.1000 

2.41145 

0.1000 

2.17337 

0.1000 

1.96925 

0.1000 

1.82852 

0.1000 

1.77779 

0.1000 

2.60044 

2.41991 

2.21757 

2.00791 

1.81024 

1.64660 

1.53806 

1.50000 
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Table  4.  Values  of  semi-major  axis  for  the  general  case,  T  =  0.6 


0  1  2  3  21  22  23  T 


15.0  0.41785  0.79604  0.83617  ... 

30.0  0.42136  0.80879  0.85630  ... 

45.0  0.42702  0.82888  0.88714  ... 

60.0  0.43458  0.85472  0.92534  ... 

75.0  0.44367  0.88438  0.96729  ... 

90.0  0.45383  0.91576  1.00966  ... 

105.0  0.46449  0.94677  1.04973  ... 

120.0  0.47500  0.97547  1.08538  ... 

135.0  0.48463  1.00011  1.11500  ... 

150.0  0.49259  1.01919  1.13732  ... 

165.0  0.49800  1.03142  1.15134  ... 

180.0  0.50000  1.03571  1.15618  ... 

195.0  0.49772  1.03114  1.15113  ... 

210.0  0.49038  1.01695  1.13564  ... 

225.0  0.47736  0.99253  1.10941  ... 

240.0  0.45833  0.95752  1.07241  ... 

255.0  0.43341  0.91186  1.02497  ... 

270.0  0.40332  0.85602  0.96778  ... 

285.0  0.36960  0.79127  0.90183  ... 

300.0  0.33465  0.72029  0.82825  ... 

315.0  0.30171  0.64797  0.74886  ... 

330.0  0.27443  0.58248  0.66912  ... 

345.0  0.25634  0.53528  0.60417  ... 

360.0  0.25000  0.51786  0.57809  ... 


0.82821  0.82821  0.82821  0.6000 
0.84994  0.84994  0.84994  0.6000 
0.88229  0.88229  0.88229  0.6000 
0.92125  0.92125  0.92125  0.6000 
0.96315  0.96315  0.96315  0.6000 
1.00498  1.00498  1.00498  0.6000 
1.04438  1.04438  1.04438  0.6000 
1.07946  1.07946  1.07946  0.6000 
1.10868  1.10868  1.10868  0.6000 
1.13078  1.13078  1.13078  0.6000 
1.14468  1.14468  1.14468  0.6000 
1.14948  1.14948  1.14948  0.6000 
1.14446  1.14446  1.14446  0.6000 
1.12903  1.12903  1.12903  0.6000 
1.10285  1.10285  1.10285  0.6000 
1.06583  1.06583  1.06583  0.6000 
1.01829  1.01829  1.01829  0.6000 
0.96103  0.96103  0.96103  0.6000 
0.89548  0.89548  0.89548  0.6000 
0.82380  0.82380  0.82380  0.6000 
0.74893  0.74893  0.74893  0.6000 
0.67482  0.67482  0.67482  0.6000 
0.60857  0.60857  0.60857  0.6000 
0.57474  0.57474  0.57474  0.6000 
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Table  5.  Values  of  semi-major  axis  for  the  general  case,  T  =  1.1 


0 

1 

2 

3 

21 

22 

23 

T 

15.0 

0.22792 

0.60611 

0.67968 

0.65760 

0.65760 

0.65760 

1.1000 

30.0 

0.22983 

0.61726 

0.70437 

0.68722 

0.68722 

0.68722 

1.1000 

45.0 

0.23292 

0.63478 

0.74160 

0.72855 

0.72855 

0.72855 

1.1000 

60.0 

0.23704 

0.65719 

0.78666 

0.77515 

0.77515 

0.77515 

1.1000 

75.0 

0.24200 

0.68272 

0.83471 

0.82246 

0.82246 

0.82246 

1.1000 

90.0 

0.24754 

0.70948 

0.88163 

0.86747 

0.86748 

0.86746 

1.1001 

105.0 

0.25336 

0.73564 

0.92440 

0.90820 

0.90812 

0.90822 

1.0994 

120.0 

0.25909 

0.75956 

0.96106 

0.94302 

0.94325 

0.94306 

1.1010 

135.0 

0.26435 

0.77983 

0.99044 

0.97144 

0.97134 

0.97137 

1.1000 

150.0 

0.26868 

0.79529 

1.01186 

0.99217 

0.99205 

0.99218 

1.0992 

165.0 

0.27164 

0.80506 

1.02491 

1.00483 

1.00480 

1.00484 

1.0997 

180.0 

0.27273 

0.80844 

1.02929 

1.00913 

1.00909 

1.00913 

1.0997 

195.0 

0.27148 

0.80491 

1.02488 

1.00478 

1.00475 

1.00478 

1.0998 

210.0 

0.26748 

0.79405 

1.01166 

0.99169 

0.99175 

0.99169 

1.1004 

225.0 

0.26038 

0.77555 

0.98983 

0.97007 

0.97000 

0.97016 

1.0986 

240.0 

0.25000 

0.74919 

0.95982 

0.94033 

0.93992 

0.94000 

1.1007 

255.0 

0.23640 

0.71486 

0.92222 

0.90156 

0.90239 

0.90251 

1.0945 

270.0 

0.21999 

0.67269 

0.87759 

0.85796 

0.85796 

0.85534 

1.1485 

285.0 

0.20160 

0.62327 

0.82596 

0.80939 

0.80155 

0.80800 

1.0538 

360.0 

0.13636 

0.40422 

0.51465 

0.50456 

0.50454 

0.50456 

1.0997 

Some  transfer  angle  cases  were  omitted  because  T  =  1.1  exceeded  the  minimum 
energy  transfer  time,  tme. 
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Table  6.  Comparison  of  Gauss  Method  with  series  solution,  T  =  0.1 
Series  Solution:  (repeated  for  convenience) 


9 

1 

2 

3 

21 

22 

23 

T 

15.0 

2.50711 

2.88530 

2.89199 

...  2.89171 

2.89171 

2.89171 

0.1000 

30.0 

2.52813 

2.91557 

2.92349 

...  2.92326 

2.92326 

2.92326 

0.1000 

45.0 

2.56211 

2.96397 

2.97368 

...  2.97351 

2.97351 

2.97351 

0.1000 

60.0 

2.60747 

3.02762 

3.03939 

...  3.03924 

3.03924 

3.03924 

0.1000 

75.0 

2.66201 

3.10273 

3.11655 

...  3.11641 

3.11641 

3.11641 

0.1000 

Gauss  Method: 


9 

1 

2 

3 

21 

22 

23 

T 

15.0 

2.67394 

2.89500 

2.89166 

...  2.89171 

2.89171 

2.89171 

0.1000 

30.0 

2.16686 

2.98185 

2.92002 

...  2.92326 

2.92326 

2.92326 

0.1000 

45.0 

1.57486 

3.35599 

2.92979 

...  2.97351 

2.97351 

2.97351 

0.1000 

60.0 

1.04310 

5.36176 

3.03939 

...  3.03924 

3.03924 

3.03924 

0.1000 

The  Gauss  method  failed  to  converge  for  6  >  75°. 


Particular  Cases 


In  order  to  numerically  integrate  the  equations  of  motion  for  an  n-body  orbit  problem, 
an  initial  estimate  of  a  trajectory  must  be  provided.  Generally,  the  two  body  solution  is 
used  to  obtain  this  trajectory  information  when  one  gravitational  source  is  predominant 
over  the  others.  The  two  body  series  solution  presented  in  the  analytical  development 
section  was  used  to  provide  the  initial  trajectories  for  two  multiple  body  problems.  The 
first  case  resulted  in  a  hyperbolic  transfer,  and  the  second  produced  an  elliptic  transfer. 

The  following  values  may  be  considered  to  be  given  data  for  the  two-body  problem. 
The  position  coordinates  are  given  in  Astronomical  Units,  {AU)  where  one  AU  is  the 
mean  distance  from  the  Sun  to  the  Earth.  The  time  unit  used  is  the  day,  ( d )  where  one 
day  equals  24  hours.  Angular  values  axe  presented  in  degrees.  The  gravitational  constant 
k  =  y/ji.  The  first  case  was  as  follows: 

xl  =  0.46918988885509  AU 

yl  =  -0.77383205171227  AU 

zl  =  -0.01964834734771  AU 

x2  =  1.31776281141600  AU 

y2  =  -0.41736193703330  AU 

z2  =  0.02991885008669  AU 

k  =  0.01720209895000  AU^/d 
t  =  40.000  d 

rl  =  0.90517470889026  AU 

r2  =  1.38260079242914  AU 

0  =  41.26824° 

T  =  -0.0059691678669 
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From  the  given  data,  the  B{  coefficients  of  equation  (23)  may  be  calculated.  A  partial 
sum  is  formed  of  the  first  »  terms,  then  the  result  is  multiplied  by  the  quantity  (s/2),  where  s 
is  the  semi-perimeter.  This  product  gives  the  semi-major  axis.  The  following  output  shows 
the  value  of  semi-major  axis  obtained  by  using  the  first  i  terms  in  the  series.  The  first 
number  is  the  index  i,  and  the  second  number  is  the  semi-major  axis  in  ( AU ). 

1  -49.2301806515044904 

2  -48.7672453986109379 

3  -48.7679313499686148 

4  -48.7679320992466774 

5  -48.7679321023198326 

6  -48.7679321023313657 

7  -48.7679321023314029 

8  -48.7679321023314030 

9  -48.7679321023314030 

Once  the  semi-major  axis  has  been  found,  the  orbit  is  completely  determined  since 
two  positions  on  the  orbit  are  known  from  the  given  data.  Any  of  the  remaining  classical 
orbital  elements  may  be  calculated  as  well  as  position  and  velocity  vectors  for  any  point 
in  the  orbit.  In  the  multibody  problems  it  was  desired  to  know  the  components  of  the 
velocity  vector  at  position  rj,  so  they  are  provided  for  reference.  These  components  were 
used  as  the  initial  value  in  a  boundary  value  problem. 

x  =  2.5256692268401 1£  -  0002  AU/d 
y  =  5.02133832285282£  —  0003  AU/d 
z  =  1.20446318655231.E  -  0003  AU/d 
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Following  the  same  format  as  before,  the  second  case  was  as  follows: 


*1  =  0.50186422427732  AU 

yl  =  -0.77640603245208  AU 

zl  =  -0.01549685878577  AU 

x2  =  1.37003894998300  AU 

y2  =  -0.21022615184980  AU 

z2  =  0.02453126302031  AU 

k  =  0.01720209895000  AU^/d 
t  =  54.000  d 

rl  =  0.92461569285281  AU 

r2  =  1.38629129055097  AU 

0  =  48.43571° 

T  =  0.1886166547276 
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i,  semi-major  axis  (AU): 


1  1.58419934014415077 

2  2.05984061297696060 

3  2.08346900567006944 

4  2.08277664547186128 

5  2.08286558797827610 

6  2.08285434745162601 

7  2.08285556050262464 

8  2.08285544531016136 

9  2.08285545548735275 

10  2.08285545459028151 

11  2.08285545467360664 

12  2.08285545466545692 

13  2.08285545466626006 

14  2.08285545466618325 

15  2.08285545466619034 

16  2.08285545466618969 

17  2.08285545466618975 

18  2.08285545466618975 

19  2.08285545466618975 

20  2.08285545466618975 

initial  velocity  components: 


x  =  2. 14961598862402.fi1  -  0002  AU /d 
y  =  5.951346004451281?  -  0003  AU jd 
z  =  7.08698265474608£  -  0004  AU/d 
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IV.  Conclusions  and  Recommendations 


Summary  and  Conclusions 

The  two  body  series  solution  presented  in  the  analytical  development  section  appears 
to  be  well  behaved  within  the  following  limits:  -1  <  T  <  Tme,  3  <  i  <  23.  Note  when 
T  =  0,  the  orbit  is  parabolic,  and  the  semi-major  axis  is  infinite.  However,  the  parabolic 
orbit  is  uniquely  defined  for  a  given  problem  goemetry. 

The  series  solution  also  appears  to  be  relatively  insensitive  to  the  transfer  angle,  6. 
Presumably  this  may  be  explained  by  noting  Lambert’s  Time  Function,  from  which  the 
series  was  derived,  is  also  relatively  insensitive  to  the  transfer  angle. 

The  series  may  be  used  for  any  two  body  orbit  flight  time  from  hyperbolic  transfers 
up  to  the  minimum  energy  elliptic  transfer.  The  series  should  not  be  used  for  flight  times 
exceeding  the  minimum  energy  transfer  time,  because  the  form  of  Lambert’s  Time  Function 
used  in  the  derivation  does  not  apply  to  such  orbits.  For  flight  times  greater  than  tme, 
Lambert’s  Time  Function  includes  a  constant  term  related  to  the  orbital  period,  which 
is  the  time  required  ^  uake  one  complete  orbital  revolution.  The  presence  of  this  term 
makes  it  difficult  tr  express  Lambert’s  Time  Function  as  a  power  series. 

Because  the  series  may  be  used  for  both  hyperbolic  and  elliptic  orbits,  it  avoids  the 
necessity  of  programming  for  separate  cases  as  is  usually  done  in  the  direct  application  of 
Lambert’s  Time  Function. 

The  greatest  advantage  gained  by  using  the  series  solution  for  semi-major  axis  over 
other  methods  is  that  no  initial  value  is  required,  because  no  root-finding  technique  is 
necessary.  Given  the  flight  time  and  orbital  geometry,  it  is  possible  to  substitute  the 
appropriate  values  in  to  the  series  definition,  then  evaluate  the  series  to  obtain  the  semi¬ 
major  axis  of  the  orbit. 
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Recommendation  for  Further  Study 

The  form  of  Lambert’s  Time  Function  for  orbital  transfers  greater  than  the  minimum 
energy  transfer  time  does  not  lend  itself  to  series  representation.  If  a  series  representa¬ 
tion  were  to  be  found,  then  the  inversion/reversion  techniques  presented  in  the  analytical 
development  section  may  prove  useful  in  solving  such  a  series. 

It  may  be  possible  to  use  Asymptotic  Matching  (5:  p.  270-279),  to  find  the  series  that 
describes  greater  that  minimum  energy  orbits.  By  allowing  the  value  of  semi-major  axis 
to  approach  infinity  in  Lambert’s  Time  Function,  one  may  obtain  a  relationship  between 
flight  time  and  semi-major  axis  for  very  long  transfer  times.  This  may  be  matched  to 
the  solution  for  transfer  times  approaching  the  minimum  energy  time.  If  the  resulting 
expression  could  be  expressed  as  a  series,  the  reversion/inversion  techniques  again  may  be 
useful  in  solving  the  series. 
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Appendix  A.  Example  Coefficients 


The  first  six  coefficients  are  given: 


=  A\ 

B2  =  A2/A\ 

B3  =  (A1A3  -  A\)/A\ 

B4  =  (A21A4-3AlA2A3  +  2A32)/Asl 
Bs  =  (A3lAs-4AlA2A4-2AlAl  +  10A1AlA3-5A42)/Al 
Be  =  (A?A6  -  5 A\A2Ae  4-  15A?A^A4  -  5 A3xA3A4  +  15A?A2Al 
— 35Ai  AjA3  +  14  A§)/A? 

Br  =  ( A\A7  -  6 A\A2AG  +  21  A\A\Ah  -  6A? A3A5  -  3 A\A\  +  42A?A2A3A4 
-56 A\A32A4  +  7A?A^  -  84A1A2A3  +  126AiA^A3  -  42 A^/A}1 

It  is  not  recommended  that  these  equations  be  programmed  directly  due  to  their 
complexity.  To  avoid  coding  errors,  it  would  be  more  efficient  to  program  the  recursive 
relationships  developed  in  the  analytical  development  section.  The  recursive  algorithm 
allows  the  calculation  of  terms  of  much  higher  order  than  the  first  seven  without  explicit 
coding.  A  sample  FORTRAN  listing  is  included  in  Appendix  C  that  includes  the  necessary 
algorithms. 
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Appendix  B.  Accuracy  Plots  for  Various  Transfer  Angles 


Plots  of  accuracy  vs.  T  and  number  of  terms  are  given  for  the  following  values 
of  transfer  angle:  6  =  90°,  180°,  270°,  360°.  The  vertical  axis  represents  the  number  of 
significant  figures  of  accuracy  in  the  calculated  value  of  semi-major  axis.  In  general,  the 
number  of  significant  figures  increases  with  the  number  of  terms  in  the  series.  The  increase 
is  largest  in  the  neighborhood  of  the  Parabolic  Transfer  Time  which  corresponds  to  T  =  0. 

An  example  point  is  shown  on  figure  3,  where  the  indicated  point  has  coordinates 
T  =  -3,  i  =  8,  and  significant  figures  =  7. 


SIG.  FIG. 
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Figure  4.  Accuracy  plot  for  180  degree  transfer  angle. 
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SIC.  FIG. 


Figure  5.  Accuracy  plot  for  270  degree  transfer  angle. 
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Figure  6.  Accuracy  plot  for  360  degree  transfer  angle. 
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Appendix  C.  FORTRAN  listing  of  Algorithm 


This  is  the  FORTRAN  listing  of  the  algorithm.  The  program  requires  two  position 
vectors,  a  gravitational  constant,  and  a  flight  time,  and  will  produce  the  corresponding 
value  of  semi-major  axis. 
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PROGRAM  ELAH9 

C 

C  THIS  VERSION  HILL  HANDLE  HYPERBOLIC  AND  ELLIPTIC  ORBITS 
C 

DOUBLE  PRECISION  Q(78,78>  ,  A(78)  ,  B(78) 

DOUBLE  PRECISION  X,  XB,  S 

DOUBLE  PRECISION  TOF  ,  TIME  ,  PARATIME  ,  SEMIMAJORAXIS 
DOUBLE  PRECISION  XI  ,  Y1  ,  21  ,  X2  ,  Y2  ,  Z2  ,  K 
DOUBLE  PRECISION  R1  ,  R2  ,  DOT  ,  ARG  ,  BETA 
DOUBLE  PRECISION  CHORD  ,  SEMI PERI METER 
DOUBLE  PRECISION  K1  ,  K2  ,  K3 

DOUBLE  PRECISION  ARG1  ,  ARG2  ,  Cl  ,  C2  ,  C3  ,  TC 
DOUBLE  PRECISION  RV1  ,  RV2  ,  RV3  ,  AL  ,  BE  ,  EA  ,  F  ,  G 
DOUBLE  PRECISION  VI  ,  V2  ,  V3  ,  II 

INTEGER  I  ,  J  ,  M  ,  N  ,  FLAG  ,  FLAG2 

INTEGER  ORDER  ,  AA  ,  BB  ,  P 

COMMON  /ELAHCOM/  Q,  A,  B,  X,  XB,  S,  TOF,  TIME,  PARATIME, 

+  SEMIMAJORAXIS,  XI,  Yl,  21,  X2,  Y2,  22,  K, 

+  Rl,  R2,  DOT,  ARG,  BETA,  CHORD,  SEMIPERIMETER, 

+  Kl,  K2,  K3,  ARG1,  ARG2,  Cl,  C2,  C3,  TC,  RV1, 

+  RV2,  RV3,  AL,  BE,  EA,  F,  G,  VI,  V2,  V3,  II, 

+  I,  J,  M,  N,  FLAG,  FLAG2,  ORDER,  AA,  BB,  P 

DOUBLE  PRECISION  PI 

PARAMETER  (PI  -  3 . 141592653589D0) 

C  DEFINE  THE  FUNCTIONS  FACT  AND  POCH 

DOUBLE  PRECISION  FACT,  POCH,  TEMP 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  BEGINNING  OF  MAIN  PROGRAM 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  GEOMETRY 


xl 

m 

0. 5018 6422 42 77 32  DO 

yi 

m 

-0.77640603245208D0 

zl 

- 

-0.01549685878577D0 

x2 

• 

1 . 37003894998300D0 

y2 

a 

-0. 21022615184980D0 

z2 

* 

0.024531 26302031 DO 

c 

xl 

• 

0. 46918988885509D0 

c 

yi 

■ 

-0. 7738320517 12 27D0 

c 

zl 

* 

-0.01 964834734771 DO 

c 

x2 

• 

1.31776281141600D0 

c 

y2 

— 

-0. 41736193703330D0 

c 

z2 

- 

0. 02991885008669D0 

C  SOLAR  GRAVITATIONAL  CONSTANT  (  AU*3/2  /  DAY  ) 

K  •  0.01720209895000D0 


Figure  7.  FORTRAN  listing  of  algorithm 
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R1  -  DSQRT (  XI  *  XI  +  Y1  *  Y1  +  21  *  21  ) 

R 2  -  DSQRT (  X2  *  X2  +  Y2  *  Y2  +  22  *  Z2  ) 

DOT  -  (XI  *  X2)  +  (Y1  *  Y2)  +  (Z1  *  22) 

BETA  -  DACOS (  DOT  /  (R1  *  R2)  ) 

IF  (BETA  .LT.  PI)  THEN 
FLAG  -  -1 
ELSE 

FLAG  -  1 
END  IF 

CHORD  -  DSQRT (R1  *  R1  +  R2  *  R2  -  2.0D0  *  R1  *  R2  *  DCOS(BETA)) 
SEMIPERIMETER  -  (R1  +  R2  +  CHORD)  /  2.0D0 

PARATIME  -  (DSQRT (2 . 0D0)  /  K)  *  (SEMIPERIMETER* *1 . 5D0  + 

+  DBLE (FLAG)  *  (SEMI PERI METER  -  CHORD) **1 . 5D0)  /  3.0D0 

K1  -  DSQRT (SEMIPERIMETER* *3 . 0D0  I  2.0D0)  /  K 

K2  -  DBLE (FLAG)  *  DSQRT ( (SEMIPERIMETER-CHORD) **3. 0D0  /  2.0DO)  /K 
K3  -  (SEMIPERIMETER  -  CHORD)  /  SEMIPERIMETER 

C  INITIALIZE  MATRIX 

ORDER  -  23 

DO  20  I  -  1, ORDER 
DO  10  J  -  1, ORDER 
Q(I,J)  -  0.0D0 


10 

CONTINUE 

A(I)  -  0 

B (I)  -  0 

20 

CONTINUE 

C  CALCULATE  COEFFICIENTS  A  I  OF  ORIGINAL  SERIES 


DO  30  I  -  1, ORDER 
TEMP  -  DBLE (I) 

A(I)  -  POCH (1 . 5D0, TEMP)  *  POCH (0.5D0rTEMP>  /  POCH (2 . 5D0, TEMP) 
A(I)  -  A(I)  *  ( (1.0D0  +  DBLE (FLAG)  *  K3**(TEMP  +  1.5D0))  / 

+  (1.0D0  +  DBLE (FLAG)  *  K3**1.5D0))  /  FACT (TEMP) 

30  CONTINUE 

C  FORM  Q  MATRIX 

Q(l,l)  -  l.ODO  /  A(l) 

DO  70  I  -  2, ORDER 
DO  50  J  -  2, ORDER 
IF  (I  .(X.  J)  THEN 
DO  40  M  -  1, (I  -  1) 

Qd.J)  -  Q(I,J)  +  Q  (I-M,  1)  *  Q(M,  (J-l) ) 

40  CONTINUE 

END  IF 

50  CONTINUE 

DO  60  N  -  2,1 

Q(I.l)  -  Qd»  1)  +  A (N)  *  Q (I, N) 

60  CONTINUE 


Figure  8.  FORTRAN  listing  of  algorithm  (continued) 


35 


70 


Q(I.1>  -  Q(I,1)  /  (-A(l) ) 

CONTI HOT 

C  OSS  A  AND  Q  TO  GET  FINAL  COEFFICIENTS 

B (1)  -  A(l) 

DO  90  I  -  2, ORDER 
DO  80  M  -  2,1 

B (I)  -  B(I)  +  A (M)  *  Q( (I-I) , (M-I) ) 
80  CONTINUE 

90  CONTINUE 

C  USE  COEFFICIENTS  WITH  DIFFERENT  TIMES 


C  INITIALIZE 

X  -  0.0D0 
XB  -  O.ODO 
S  -  O.ODO 

TOF  -  S4.0D0 

TIME  -  TOF  /  PARATIME  -  1.0D0 

DO  300  I  -  1, ORDER 

XB  -  B (I)  *  TIME** (1-2) 

X  -  X  +  XB 

S  -  X  *  SEMIPERIMETER  /  2.0D0 
PRINT  *,  I,  S 

300  CONTINUE 

STOP 

END 

DOUBLE  PRECISION  FUNCTION  FACT(U) 

C  CALCULATES  FACTORIALS 

DOUBLE  PRECISION  U,  PROD,  TEMP 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

C  BEGINNING  OF  FUNCTION  FACT 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

IF  (0  .EQ.  O.ODO)  THEN 
FACT  -  1.0D0 
ELSE 

PROO  -  U 

TEW  -  U 

10  CONTINUE 

IF  (TEMP  .GT.  1.0D0)  THEN 
TEMP  -  TEW  -  1.0D0 
PROD  -  PROO  *  TEMP 


Figure  9.  FORTRAN  listing  of  algorithm  (continued) 
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onn 


GOTO  10 

caoir 

FACT  -  PROD 
ENDir 

END 

DOUBLE  PRECISION  FUNCTION  P0CH(U,V> 

C  CALCULATES  POCHHAMMER  SYMBOLS 

DOUBLE  PRECISION  U,V, PROD, TEMP 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

BEGINNING  OF  FUNCTION  POCH 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

PROD  -  1.0D0 
TEMP  -  V 

10  CONTINUE 

IF  (TEMP  .GE.  1.0D0)  THEN 
TEMP  -  TEMP  -  1.0D0 
PROD  -  PROD  *  (U  +  TEMP) 

GOTO  10 
END  IF 

POCH  -  PROD 
END 


Figure  10.  FORTRAN  listing  of  algorithm  (continued) 
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